function fig = plotTfcPeaks_Core(datatable, header, subjectname, epochDatasetName, cfg)

%% Start
disp('Plotting figures...')
if cfg.PlotVisible
    fig = figure();
else
    fig = figure('Visible', 'off');
end
hold on;

%% Filter the data
tfSelect = true(size(datatable,1),1);
% Select the significant peaks
if ~isnan(cfg.PValueLessThanSelection)
    tfSelect = tfSelect & ([datatable{:,strcmp(header,'P')}]<cfg.PValueLessThanSelection)';
end
%Find the electrodes marked as selected in config.
if ~isempty(cfg.ElectrodeSelection)
    tfChSel = false(size(datatable,1),1);
    for i = 1:length(cfg.ElectrodeSelection)
        tfChSel = tfChSel | strcmp(cfg.ElectrodeSelection{i},{datatable{:,strcmp(header,'Channel')}});
    end
    tfSelect = tfSelect & tfChSel;
end
% Select the correct condition
if ~isempty(cfg.ConditionSelection)
    tfSelect = tfSelect & strcmp({datatable{:,strcmp(header,'Condition')}},cfg.ConditionSelection)';
end
% Execute the filter
datatable = datatable(tfSelect,:);

%% Preprocess the plot
% Get the list of channels from the datatable
electrodes = unique({datatable{:,strcmp(header,'Channel')}}');

% Find the regions and sort them according to the sequence in cfg.
regions = unique({datatable{:,strcmp(header,'Region')}}');
selectionOfRegions = {};
if ~isempty(cfg.RegionSelection)
    for i = 1:length(regions)
        if any(strcmp(cfg.RegionSelection, regions{i}))
            selectionOfRegions{end+1} = regions{i};
        end
    end
    regions = selectionOfRegions;
end
regionsSort = 999*ones(size(regions));
for i = 1:length(cfg.PreferredOrdering.Regions)
    Idx = strcmp(regions, cfg.PreferredOrdering.Regions(i));
    regionsSort(Idx) = i;
end
[XXX,I] = sort(regionsSort);
regions = regions(I); % Translate the regions to the preferred sorted format

% Get the min/max values for the subplots
V = cell2mat(datatable(:,strcmp(header,'Value')));
minVal = floor(min(V)*5)/5;
maxVal = ceil(max(V)*5)/5;

% Find the frequency bins and set up the subplots
% [lfbB lfbIdx] = sort([datatable{:,strcmp(header,'FrequencyMin')}]); % Find the lower frequency boundaries and sort
[bins binsM] = unique({datatable{:,strcmp(header,'Bin')}}');
[lfbB lfbIdx] = sort([datatable{binsM,strcmp(header,'FrequencyMin')}]); % Find the lower frequency boundaries and sort
bins = bins(lfbIdx); % Resort the bins
binTf = false(length(bins),size(datatable,1));

for b = 1:length(bins)
    binTf(b,:) = strcmp(bins(b), {datatable{:,strcmp(header,'Bin')}});
    subplot(2,ceil(length(bins)/2),b);
    hold on;
    ylabel(bins{b});
    set(gca, ...
        'YGrid','on',...
        'XGrid','on',...
        'XLim',[-4000 4000],...
        'YLim',[minVal maxVal],...
        'XTick',[-4000, -3000, -2000, -1000, 0, 1000, 2000, 3000, 4000],...
        'XTickLabel',{'-4s','','-2s','','M','','2s','','4s'});
%         'XTickLabel',{'','','','','','','','',''});
    if b==2
        title([subjectname ' ' epochDatasetName]);
    end
    if b==length(bins)
    set(gca, ...
        'YGrid','on',...
        'XGrid','on',...
        'XTick',[-4000, -3000, -2000, -1000, 0, 1000, 2000, 3000, 4000],...
        'XTickLabel',{'-4s','','-2s','','M','','2s','','4s'});
    end
end

%% Add the data
LatColumn = {datatable{:,strcmp(header,'Latency (ms)')}};
ValColumn = {datatable{:,strcmp(header,'Value')}};
% lDesc = {};
mCfgMs = fieldnames(cfg.Markers);
for r = 1:length(regions)
    region = cell2mat(regions(r));
    markerIdx = [];
    for i = 1:length(mCfgMs)
        if strcmp(cfg.Markers.(mCfgMs{i}).Description, region)
            markerIdx = mCfgMs{i};
            break;
        end
        if isempty(markerIdx); markerIdx = 'Other'; end;
    end
    markerType = cfg.Markers.(markerIdx).Marker;
    markerEdgeSize = cfg.Markers.(markerIdx).MarkerEdgeSize;
    markerSize = cfg.Markers.(markerIdx).MarkerSize;
    markerEdgeColor = cfg.Markers.(markerIdx).MarkerEdgeColor;
    markerFaceColor = cfg.Markers.(markerIdx).MarkerFaceColor;
    
    for b = 1:length(bins)
        subplot(2,ceil(length(bins)/2),b);
        rtf = strcmp(region,{datatable{:,strcmp(header,'Region')}}); % find the rows for that region..
        btf = binTf(b,:); 
        X = cell2mat(LatColumn(rtf & btf));
        Y = cell2mat(ValColumn(rtf & btf));

        if isempty(markerFaceColor)
        SOpen = scatter(X,Y, ...
            'Marker', markerType, ...
            'SizeData', markerSize, ...
            'LineWidth', markerEdgeSize, ...
            'MarkerEdgeColor', markerEdgeColor);
            set(SOpen, 'DisplayName', region);
        else
        SFilled = scatter(X,Y, 'filled', ...
            'Marker', markerType, ...
            'SizeData', markerSize, ...
            'MarkerEdgeColor', markerEdgeColor, ...
            'MarkerFaceColor', markerFaceColor);
            set(SFilled, 'DisplayName', region);
        end
%         lDesc = [lDesc;[region ' (' bins{b} ')']];
    end
end
%         legend(lDesc);

subplot(2,ceil(length(bins)/2),1);
legend1 = legend('show');
set(legend1,'Position',[0.7 0.25 0.06 0.15]);

% disp('Saving figures...')
% saveFigure(fig, [getpref('ECoG','Results') 'TFC\Significant Peaks\SelectionOfElectrodes\'], [subjectinfo.PrimaryPrefix ' ' EpochDatasetName ' - SignificantPeaksByRegionPerBin - Two-sided - Selection of Electrodes']);
% close(fig);
end